function [t_out, y_out, te, ye, ie] = my_ode45(odefun, tspan, y0, options, eventfun)
    % 参数解析和初始化
    if nargin < 4
        options = struct();
    end
    if nargin < 5
        eventfun = [];
    end
    
    % 解析选项参数
    RelTol = getoption(options, 'RelTol', 1e-3);
    AbsTol = getoption(options, 'AbsTol', 1e-6);
    MaxStep = getoption(options, 'MaxStep', diff(tspan)/10);
    InitialStep = getoption(options, 'InitialStep', []);
    
    % 初始化输出变量
    t_out = tspan(1);
    y_out = y0(:)';
    te = [];
    ye = [];
    ie = [];
    
    % 事件处理初始化
    have_events = ~isempty(eventfun);
    if have_events
        [v, isterm, dir] = eventfun(tspan(1), y0);
        n_events = length(v);
        v_old = v;
        isterm = isterm(:);
        dir = dir(:);
    else
        n_events = 0;
        v_old = [];
        isterm = [];
        dir = [];
    end
    
    % 初始化时间和状态
    t = tspan(1);
    y = y0(:);
    h = compute_initial_step(odefun, t, y, tspan(2), RelTol, AbsTol, MaxStep, InitialStep);
    
    % 主循环
    while t < tspan(end)
        h = min(h, tspan(end) - t);
        
        % 计算一步RK
        [k, y_new, y_hat, err] = rk_step(odefun, t, y, h);
        
        % 计算误差
        err = max(abs(err) ./ (AbsTol + RelTol * max(abs(y), abs(y_new))));
        if err <= 1
            % 接受步长
            t_new = t + h;
            
            % 处理事件
            if have_events
                [v_new, ~, ~] = eventfun(t_new, y_new(:)');
                events = check_events(v_old, v_new, dir);
                for i = 1:n_events
                    if events(i)
                        % 插值查找事件时间
                        [tevent, yevent] = find_event(odefun, eventfun, t, y, t_new, y_new, k, h, i, v_old(i), v_new(i), dir(i), AbsTol, RelTol);
                        if ~isempty(tevent)
                            te = [te; tevent];
                            ye = [ye; yevent'];
                            ie = [ie; i];
                            % 检查是否终止
                            if isterm(i)
                                % 将事件点加入输出
                                t_out = [t_out; tevent];
                                y_out = [y_out; yevent'];
                                return;
                            end
                        end
                    end
                end
                v_old = v_new;
            end
            
            % 存储结果
            t_out = [t_out; t_new];
            y_out = [y_out; y_new(:)'];
            
            % 更新状态
            t = t_new;
            y = y_new;
            
            % 调整步长
            h = h * min(4, max(0.1, 0.8 * (1/err)^(1/5)));
        else
            % 拒绝步长，调整
            h = h * max(0.1, 0.8 * (1/err)^(1/5)));
        end
        h = min(h, MaxStep);
    end
end

function h = compute_initial_step(odefun, t, y, t_end, RelTol, AbsTol, MaxStep, InitialStep)
    if ~isempty(InitialStep)
        h = InitialStep;
        return;
    end
    h = min(t_end - t, MaxStep);
    k1 = odefun(t, y);
    if any(abs(k1) > 1e-10)
        k2 = odefun(t + h/5, y + h/5*k1);
        scale = AbsTol + RelTol * abs(y);
        err = norm(h*(k1 - k2)/5 ./ scale, inf);
        if err > 0
            h = 0.8 * h * (1 / err)^(1/5);
        end
    end
    h = min(max(h, 1e-4*(t_end - t)), MaxStep);
end

function [k, y_new, y_hat, err] = rk_step(odefun, t, y, h)
    % Dormand-Prince系数
    c = [0; 1/5; 3/10; 4/5; 8/9; 1; 1];
    A = [0,0,0,0,0,0,0;
         1/5,0,0,0,0,0,0;
         3/40,9/40,0,0,0,0,0;
         44/45,-56/15,32/9,0,0,0,0;
         19372/6561,-25360/2187,64448/6561,-212/729,0,0,0;
         9017/3168,-355/33,46732/5247,49/176,-5103/18656,0,0;
         35/384,0,500/1113,125/192,-2187/6784,11/84,0];
    b = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0];
    b_hat = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40];
    
    k = zeros(numel(y),7);
    k(:,1) = odefun(t, y);
    k(:,2) = odefun(t + c(2)*h, y + h*(A(2,1)*k(:,1)));
    k(:,3) = odefun(t + c(3)*h, y + h*(A(3,1)*k(:,1) + A(3,2)*k(:,2)));
    k(:,4) = odefun(t + c(4)*h, y + h*(A(4,1)*k(:,1) + A(4,2)*k(:,2) + A(4,3)*k(:,3)));
    k(:,5) = odefun(t + c(5)*h, y + h*(A(5,1)*k(:,1) + A(5,2)*k(:,2) + A(5,3)*k(:,3) + A(5,4)*k(:,4)));
    k(:,6) = odefun(t + c(6)*h, y + h*(A(6,1)*k(:,1) + A(6,2)*k(:,2) + A(6,3)*k(:,3) + A(6,4)*k(:,4) + A(6,5)*k(:,5)));
    k(:,7) = odefun(t + c(7)*h, y + h*(A(7,1)*k(:,1) + A(7,2)*k(:,2) + A(7,3)*k(:,3) + A(7,4)*k(:,4) + A(7,5)*k(:,5) + A(7,6)*k(:,6)));
    
    y_new = y + h * sum(bsxfun(@times, b, k), 2);
    y_hat = y + h * sum(bsxfun(@times, b_hat, k), 2);
    err = y_new - y_hat;
end

function events = check_events(v_old, v_new, dir)
    events = false(size(v_old));
    for i = 1:numel(v_old)
        if dir(i) > 0
            events(i) = (v_old(i) > 0 && v_new(i) < 0) || (v_old(i) < 0 && v_new(i) > 0);
        elseif dir(i) < 0
            events(i) = (v_old(i) < 0 && v_new(i) > 0) || (v_old(i) > 0 && v_new(i) < 0);
        else
            events(i) = (v_old(i) * v_new(i)) <= 0;
        end
    end
end

function [tevent, yevent] = find_event(odefun, eventfun, t0, y0, t1, y1, k, h, event_idx, v0, v1, dir, AbsTol, RelTol)
    tevent = [];
    yevent = [];
    if (v0 * v1) > 0
        return;
    end
    if dir > 0 && ~(v0 > 0 && v1 < 0)
        return;
    end
    if dir < 0 && ~(v0 < 0 && v1 > 0)
        return;
    end
    
    % 使用布伦特方法求解事件时间
    theta0 = 0;
    theta1 = 1;
    tol = 1e-8;
    max_iter = 100;
    for iter = 1:max_iter
        theta = (theta0 + theta1)/2;
        t_theta = t0 + theta*h;
        y_theta = interpolate(theta, k, h, y0);
        v_theta = eventfun(t_theta, y_theta');
        v_theta = v_theta(event_idx);
        
        if abs(v_theta) < tol
            break;
        end
        
        if v0 * v_theta < 0
            theta1 = theta;
            v1 = v_theta;
        else
            theta0 = theta;
            v0 = v_theta;
        end
    end
    tevent = t0 + theta*h;
    yevent = interpolate(theta, k, h, y0);
end

function y = interpolate(theta, k, h, y0)
    % Dormand-Prince插值多项式
    b = [1, 0, 0, 0, 0, 0, 0; % 线性插值示例，需替换为实际系数
         -3/2, 3, -3/2, 0, 0, 0, 0]; % 示例系数，实际需根据Dormand-Prince方法调整
    y = y0 + h * sum(bsxfun(@times, b(theta+1,:), k), 2);
end

function val = getoption(options, name, default)
    if isfield(options, name)
        val = options.(name);
    else
        val = default;
    end
end